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ABSTRACT 

^A  study  is  made  of  the  step-by-step  integration  of  ordinary  differential 
equations  of  the  first  order  by  means  of  formulas  obtained  from  the  Gregory- 
Newton  backward  interpolating  formula.  Tables  of  relevant  constants  are 


presented.. 
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Consider  the  ordinary  differential  equation  of  the  first  order 

y1  =  £(y,  x)  (1) 

which  is  to  beinte  grated  step-by-step  over  x  =  a  (  h  )  b  .  A  natural  method 

for  accomplishing  this  integration  is  a  predictor-corrector  process  based  upon 

a  suitable  finite  difference  interpolating  formula. 

Let  x  be  a  tabular  point  within  (a,  b)  and  assume  that  y  0is  known 
o  7 

x  -  a 

at  the  points  x  j  s  xq  -  jh  ( j  =  0,  1  ,  - - )  .  One  can  then  write  for 

y*  the  approximation  to  it  which  is  provided  by  the  Gregory -Newton  backward 
interpolating  formula 


y'=^  -jr  (u  +j  -  I)[j]  +  hj+1  -pijyr  (u+J)tJ+l]  (J+1)(?)  (2) 

j=0 

where  u  =  (x  -  xo)/h,  f  j  is  the  forward  difference  of  yr  about 

x  j  >  x  J  —  §  <  xq,  and  (u  -  1/)^  is  the  factorial  polynomial 

(  u-i)UL(u-i)  (u -i-1)-  •  •  (u-i-j  +  1)  (3) 

which  possesses  the  expansion 


(u.  1)^=2  iskuj_k  m 

k=0 

where  the  are  the  generalized  Stirling  numbers  of  the  first  kind  [l]  . 

The  formula  (2)  can  be  integrated  in  two  fashions.  In  the  first  it  is  assumed 
that  y^  ,  y  ^  ,  .  .  .  ,  y_  j  are  accurately  known  and  that  a  prediction  of 
y^  =  y(x^  )  is  desired: 
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J 

V 


y(x1)  =  yo  +  0. 

j=0 


PREDICTOR 


(4) 


where  the  error  Ej  is  given  by 


and 


■Ej|  <  hJ+2  pj+1  |f(J+1)| 


max 


Pj  =  -yr  j  (u+j  -Dtj]  du 


(5) 


(6a) 


ci^ 


j  Zu  j  + 1  -P 

p=0 

L  ( j )  j  • 


-(j-D 


(6b) 


(6c) 


where  L(j)  is  the  least  common  denominator  of  l/l,  1  /Z,  .  .  .  ,  l/(j  +  1). 

In  the  second  method  it  is  assumed  that  y  p  .  .  .  ,  y  j  are  accurately  known, 
that  y^  is  approximately  known,  and  that  a  corrected  value  of  yQ  is  desired: 

J 

y(x  )  =  y  .  +  )  p*  A^'f  .  +  E*  CORRECTOR  (7) 

°  “  1  l—L  J  •  -J  J 

j=0 

* 

where  the  error  Ej  is  given  by 

,|E;i<  hJ+2i,3J+1*iif<j+i* 


max 


(8) 
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and 


0 

P*  =  IT  J  (u+j-  D[j]  du 
-  1 


(9a) 


-  3 

"  "LTjTjT 


j  + 1  -p 


s  j 

p 


(9b) 


(9c) 


This  corrected  value  of  can  itself  be  used  in  (7)  to  obtain  what  is 
presumably  a  still  better  value  of  y^  ,  and  this  process  can  be  repeated 
indefinitely  until  it  converges  to  a  final  value  or,  in  bad  cases,  is  seen  to  be 
divergent.  The  use  of  (4)  to  obtain  an  initial  estimate  and  the  repeated 
use  of  (7)  to  improve  this  value  constitutes  a  well-known  predictor- corrector 
process , 

The  quantities  j5j  and  fk  ,  as  defined  by  (6a)  and  (9a)  ,  respectively, 
have  been  described  by  Collatz  [2]  and  tabulated  by  him  for  j  =  0(1)6.  The 
increasing  use  of  digital  machines  -  often  in  double  precision  -  for  the 
integration  of  differential  equations  has  made  a  somewhat  extended  table 
of  these  coefficients  desirable  and  the  existence  now  [l]  of  extensive  tables 
of  the  has  made  possible  the  simple  computation  from  (6b)  and  (9b). 

Table  I  presents  the  results  of  such  an  extension;  the  table  entries  were 
checked  using  the  recursion  relation  [2] 


Pj+i  Pj +  Pj+i 


(10) 
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Figure  I  presents  a  graph  of  these  coefficients.  The  p.  can  be  seen  to  fall 

J 

off  slowly  and  can  be  shown,  from  (6a),  to  satisfy  the  inequality 


<  1 


j>0  (11a) 


which  implies  that  their  decline  is  very  slow  indeed.  The  p^  can  be  seen 
to  fall  off  somewhat  faster;  this  is  to  be  expected  since,  from  (6a)  and  (9a), 


A 


j  >2  (11b) 


These  data  point  up  the  intrinsic  superiorities  of  corrector  formulas  over 
predictor  formulas :  (i); that,  since  |  p^*  |  <  I P j  I  the  series  (7)  converges  faster 
than  the  series  (4),  and  (i  i  )  that,  since  the  p^  decay  much  faster  than  the 
p.  ,  the  buildup  of  computational  error  in  the  taking  of  successive  differences 

J  * 

will,  for  a  given  number  of  terms  in  the  interpolating  series,  have  much 
less  effect  on  the  final  answer  when  a  corrector  formula  is  used. 

In  actual  practice  the  calculation  of  the  several  differences  is  often 
not  carried  out.  Instead, the  differences  are  expanded  as 


p=0 


and  (4)  rewritten  as 

y  (xi ) =  y0  +  h  £  £  rpj  Pj  f.p  +  Sj 

j=0  p=0 


(12) 


(13a) 


FIGURE 


* 


382  895  428  860  359  262  650  880 
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J 

V 

=  y  +  h  /  a  ( J )  f  +Et 

o  Lu  p  -p  J 

P“0 

J 

=  y»  +  lttttt-  I  aP(J)f-p  +  ej 

*  P^O 


while  (7)  can  be  rewritten  as 

J 


y(xQ)  =  y 


1 


$  # 
y  .  (3.  f  +  Et 
PJ  ^3  -P  J 


=  y 


j=o  P=o 

j 

* 

a  (  J  )  i  r  T 

p  -p  J 


-i + » i  v<,,f-+  e-* 


(13b) 


(13c) 


(14a) 


(14b’ 


p=0 


-  M  +  -LTJPT-  X  5P*  +  EJ*  .  <lte» 

P~0 

The  calculation  of  the  next  value  of  y  can  then  be  accomplished  directly 

from  (13c)  or  (14c)  and  the  labor  of  maintaining  a  difference  table  thereby 

* 

eliminated.  Tables  of  6^  (  J)  and  6^  (  J)  have  been  computed  for  J  =  0(1)10 
and  p  =  0(1  )J  and  are  given  in  Tables  II  and  III?  respectively*  they 


represent  a  considerable  extension  over  the  existing  tables  [2,  3]  which  go 
at  most  to  J  =  5 
the  relations  [2] 


at  most  to  J  =  5  .  The  values  of  6  (  J)  and  6  (  J)  were  checked  by 

p  p 


The  values  of  a  (5)  given  in  Reference  3  are  believed  to  be  in  error. 
P 
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(15a) 


(15b) 


and  by  having  key  portions  of  the  computations  individually  repeated, 

The  selection  of  a  predictor-corrector  pair  for  a  specific  problem 
is  by  no  means  simple,  it  being  necessary  to  choose  J  and  h  with  care 
to  minimize  the  effects  of  roundoff  and  truncation  error.  However,  for 
modern  digital  machines  operated  in  double  precision,  J  =  10  will  probably 

be  as  large  as  is  profitable  since  by  this  point  the  computing  error  generated 

i  # 

in  calculating  the  A  f  .  will  be  growing  much  faster  than  the  (3 .  will  be 

J  J 

decreasing,  and  the  total  error  due  to  this  cause  will  normally  be  at  least  as 
important  as  the  truncation  error  Ej.  ,  Finally  ,  for  suitable  f(y,  x),  it 
may  be  desirable  to  suspend  the  common  practice  of  using  a  low  order  predictor 
to  obtain  an  estimate  of  the  next  point  which  can  then  be  improved  using  a 
high  order  corrector,*  the  existence  of  high  order  predictors  makes  the 
initial  use  of  a  high  order  predictor  formula  and  the  suspension  of  corrector 
iterations  seem  attractive  where  the  superior  roundoff  suppressing  properties 
of  a  corrector  formula  are  not  essential. 


c? 
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